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ABSTRACT 

General analytic arguments lead us to expect that in the modified dynamics 
(MOND) self-gravitating disks are more stable than their like in Newtonian 
dynamics. We study this question numerically, using a particle-mesh code 
based on a multi-grid solver for the (nonlinear) MOND field equation. We 
start with equilibrium distribution functions for MOND disk models having 
a smoothly truncated, exponential surface- density profiles and a constant 
Toomre Q parameter. We find that, indeed, disks of a given "temperature" are 
locally more stable in MOND than in Newtonian dynamics. As regards global 
instability to bar formation, we find that as the mean acceleration in the disk 
is lowered, the stability of the disk is increased as we cross from the Newtonian 
to the MOND regime. The degree of stability levels off deep in the MOND 
regime, as expected from scaling laws in MOND. For the disk model we use, 
this maximum degree of stability is similar to the one imparted to a Newtonian 
disk by a halo three times as massive at five disk scale lengths. 



Subject headings: galaxies: kinematics and dynamics 
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1. Introduction 

Underlying MOND is the assumption that galaxies do not posses a significant dark 
halo. As pointed out by [Ustriker and Peebles (1973)| , a massive dark halo may be an 



important stabilizing agent of galactic disks. It is thus interesting to compare the stability 
of bare disks in MOND to the stability of similar Newtonian disks with dark halos. 
Such considerations may also provide a MOND explanation (see [Milgrom 1989|) of the 
revised Freeman law whereby the distribution of central surface brightnesses of galactic 
disks appears to be cutoff rather sharply above a certain surface brightness -Bq, (see e.g. 
McGaugh 1996| for a recent review and references). Translating this value of Bq into a mean 



surface density (for exponential disks) one obtains a limiting surface density that is nearly 
Eq = aoG~^, with oq the acceleration constant of MOND. In MOND, disks with a mean 
surface density S >> Sq have a different dynamical behavior than those with E < Sq. In 
particular the former are Newtonian and thus are beset by the well-known instabilities of 
bare Newtonian disks. The latter are more stable locally (as shown in [Milgrom 1989| using 



perturbation theory) and, as we shall show in this work, are also more stable globally. 
Global added stability is also supported by preliminary N-body calculations carried out by 
Christodoulou (1991)| , and by |Griv and Zhytnikov (199^ The Freeman law, which asserts 



that the former type of disk is rare, may result from this disparity. 

Toomre showed that in Newtonian dynamics a disk is stable to all local axisymmetric 
disturbances at radius R if the dimensionless quantity 

where aR is the radial velocity dispersion, k is the epicycle frequency, and S is the surface 
density. The criterion in MOND is obtained by simply replacing G by G/yU+(l + L+)^/^, 
where /i"*" is the value of the MOND interpolating function, /x, just above the disk, and 
L = din fi{x)/d ln{x) ([Milgrom 1989|) . Although Toomre's criterion rests on local analysis. 



it is found empirically that the condition Q > 1 everywhere in the disk is a necessary and 
sufficient condition for global axisymmetric stability. Stellar disks are always stable to local 
non-axisymmetric disturbances ([Goldreich and Lynden-Bell 1965[ , [Julian and Toomre 1966[ ). 



Numerical simulations have shown that stellar disks are subject to global non-axisymmetric 
instabilities, especially the bar instability. This result was confirmed analytically for a 
few models by linear, normal-mode analysis. The majority of rotationally supported, 
self-consistent disk models studied to date by numerical simulations and analytical global 
analysis are violently unstable to bar formation. However, these simulations do not reveal 
the mechanism of the instability nor suggest a way to avoid it. [Toomre (1981) suggested a 



mechanism for the bar instability based on what he named swing amplification. 
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Even the Newtonian-plus-dark-matter case the stabihty problem is anything but 
resolved. So, we shall not focus our work on testing for absolute stability in MOND. Instead 
we perform a comparative study between the added stability given to the disk by MOND, 
and that given by a dark halo. In particular we shall ask to what extent can MOND replace 
the halo's contribution to the stability of disk galaxies. 

Even for Newtonian gravity one lacks simple analytic equilibrium solutions of the 
coUisionless Boltzmann equation for a thin disk. Some of the equilibrium models studied to 
date are the Isochrone by [Kalnajs (1978)| , Kuzmin-Toomre (|Sellwood 1986| , [Hunter 1992|) , 



and the Sawamura disks ( ^awamura 1983| ). The extent to which the bar instability and 
others depend on the specific properties of these models is unknown. The situation for 
MOND is even more difficult since we have no analogous analytical models. The analytical 
methods used in [Kalnajs (1972)| and [Kalnajs (1977)[ for linear, normal-mode analysis are 



very cumbersome and give no physical insight into the nature of the instability. A simpler 
way to get the unstable modes is through N-body simulations. We have developed a 
three-dimensional, N-body code and potential solver for the nonlinear MOND problem, in 
which the potential is determined from the equation proposed by [Bekenstein and Milgrom 
(1984)1 



2. Description of the MOND potential solver 

Bekenstein and Milgrom (1984)' have formulated a non-relativistic Lagrangian theory 
for MOND, in which the acceleration field g produced by a mass distribution p is derived 
from a potential (g = — V0) satisfying the equation 

V- [/i(|V0|/ao)V0] =47rG'p (2) 

instead of the usual Poisson equation V ■ V0 = inGp, where ~ x for x <^ 1, 
and p{x) ~ 1 for x ^ 1, and ao is the acceleration constant of MOND. The form 

= x/vm? has been used in all rotation curve analyses and we also use it here. A 
solution to the field equation exists and is unique in a domain D in which p(r) is given and 
on the boundary of which (p, or yu(|V0|/ao)i9„0, is specified ( [Milgrom 1986| ). In this theory, 
the usual conservation laws of momentum, angular momentum, and energy (properly 
defined) hold, and, in addition, the center-of-mass acceleration of a star or a gas cloud in 
the field of a galaxy obeys the basic MOND assumption even if its internal accelerations are 
high. 

The nonlinearity of the MOND field eq.(^ prevents one from using the standard 
potential solvers (force calculators), at least in a straightforward way. We wrote a multigrid 
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solver for the finite difference approximation of the MOND field equation and incorporated 
it into an N-body code using the particle-mesh algorithm. We give a brief description of the 
N-body code in Appendix |B|, and also describe there some of the tests we have performed 
to establish its accuracy, and the setup of initial conditions for the simulation. We lack 
analytical potential density pairs for disks in MOND (apart for that for the Kuzmin disk), 
not to mention self-consistent stationary models. We have thus developed a numerical 
scheme for generating self-consistent stationary disk models with specified potential, surface 
density, and radial velocity dispersion. This scheme is described in appendix 0. 

Because the potential solver is novel we describe it here briefly. The discretization 
scheme used is depicted in Figure 

It uses central differencing between neighboring grid points to approximate the 
divergence and the components of V0 appearing in eq.(|]). Only for some of the components 
of V0 appearing in the argument of the function fi do we use central differencing between 
grid points that are two grid spacing apart. The part of the divergence in eq.(0) at 
point A is approximated by [S{M) — S{L)]/h , where S = fi{\V(t>\/aQ)dx(p. dx(f>{M) is 
approximated by [0(5) - (j){A)]/h, dy(f){M) by [0(/) + (j){H) - (f){J) - 0(K)]/(4/i), and 
dz4>{M) by [0(C) + (f){D) — (f){E) — 0(F)]/ (4/i). A similar calculation is done for point L and 
for the remaining parts of the divergence. This is a stable second order discretization, which, 
importantly, is flux conserving. The MOND equation, like Poisson's, can be transformed 
using Gauss's theorem into a flux equation 



fi—dS 

dD on 



AtcG / p£x, 



D 



(3) 



where D is any domain where the MOND equation is satisfied, dD is the boundary of D 
and d(f)/dn is the normal derivative of 0. The flux leaving a cell through one of its sides 
should be equal to the flux entering its neighboring cell; flux conservation means that the 








Fig. 1. — The discretization of the MOND equation. 
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two will have the same approximation in the discrete equation. 

We use the multigrid techniques developed by Brandt and collaborators ( [Brandt 



19771 , [Brandt 1984| , [Brandt 1991[ , [Bai and Brandt 1987|) , which is extremely efficient in 



solving elliptic, partial differential equations. We use the so-called fuU-multigrid algorithm 
together with the full-approximation scheme, and use Gauss-Seidel relaxation for solving the 
system of nonlinear equations produced by the discretization. Instead of solving directly for 
the new value of the unknown at the current grid point we carry out a single iteration of the 
Newton-Raphson method for finding the root of a nonlinear equation, where the derivative 
of the left-hand side of the equation with respect to a change in the unknown is calculated 
numerically. For solving the standard Poisson equation we use Gauss-Seidel relaxation with 
red-black (RB) ordering, which has two important properties: ffist, the smoothing rate for 
the usual seven-point Laplacian is the best; second, the "red" and the "black" points are 
independent and can be relaxed simultaneously. This last property is very useful in writing 
a code that is highly vectorizable and parallelizable. In order to maintain this property of 
independence in the case of the more complicated MOND equation we use a generalization 
of RB ordering using eight colors instead of two. 

The MOND potential solver was tested extensively against cases for which exact results 
are known. These include a. the complete potential field of a Kuzmin disk ( [Brada and[ 
Milgrom 1995[) ; b. the (deep) MOND, two-body force for arbitrary masses, and the N-body 



MOND force in certain symmetric configurations ([Milgrom 1994[ ); c. a general relation 
that exists between the total mass and the root-mean-square velocity for disks in the deep 
MOND regime, ffist discovered by our numerical calculations, and then proven exactly 
( [Milgrom 1994]) . 



3. Models and results 

As stated above, we concentrate on a comparative study between the stabilizing effects 
of MOND and those of dark matter halos. We have used models that have a smoothly 
truncated, exponential surface density. The disk extends out to radius Rcut = 1 (chosen as 
our unit of length), with a scale length of 0.2 in these units. The surface density is of the 
form 

E(i?) = Soexp(-i?/0.2)(l-i?^), R<Rcut (4) 

The smooth truncation of the disk is used in order to avoid the edge instabilities discussed 
by [Toomre (1981)[ , which result from a sharp drop in the surface density. We work in units 



where G = 1, = 1, and the mass is given in units of aoRl^^/G. We have constructed a 
series of models with a total mass of 0.005, 0.01, 0.02, 0.04, 0.08, 0.16, 0.32, 0.64, 1.28. The 
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disk with the lowest mass is fully in the MOND regime (a < Oq), while the disk with a 
mass of 1.28 is Newtonian almost all the way to its outer edge. The magnitude of the total 
acceleration just above the surface of the disk as a function of radius for the different mass 
models is shown in Figure ^ (This differs from the mid-plane acceleration, which enters the 
rotation curve, because of the perpendicular component that appears just above the disk.) 

A self-consistent model for a given mass distribution is also characterized by its 
"temperature": the fraction of the total kinetic energy that is in random motion. A 
convenient parameter for measuring this is the famous t = T/\W\ parameter, where 
T = / pvo^d^x is the rotational kinetic energy and \W\ is the absolute value of the 
self-gravitational energy, which by the virial relation is equal twice the total kinetic 
energy (rotational plus random) of the stationary system. In MOND, we replace the 
self-gravitational energy, in the definition of t, with twice the rotational kinetic energy of a 
cold system (where all the particles are on circular orbits). The maximum value of t is 0.5 
(realized for a cold system). The lower the value of t, the greater the part of the kinetic 
energy in random motion. Motivated by the analytical results for the Maclaurin and for 
the Kalnajs disks, and by their own numerical simulations, Ustriker and Peebles (1973] 



suggested the approximate empirical stability criterion against bar formation t < 0.14. 
Although the physics of the bar instability is only indirectly related to T/|Vr|, numerical 
studies have shown that this Ostriker-Peebles criterion provides a surprisingly useful 
empirical guide for identifying systems that are likely to be unstable. 

As a preliminary step in our comparative study we generated three self-consistent 
models for each of the total disk masses listed above. The models have a radius-independent 
value of the MOND stability parameter [see below eq.(|lD], with Q = 1.5,2.0,2.5. The 
calculated runs of t values for these models are displayed in Figure 

We then constructed MOND equilibrium models for comparative, N-body simulations. 
These were taken as smoothly truncated, exponential disks with t = T/\W\ fixed at 0.31, 
and a value of Q that is R independent. These models fall on the horizontal line in Figure 
^. While the potential field is computed on a 3-dimensional Cartesian grid, disk particles 
are, at all times, confined to the mid-plane. 

Each model was run once using MOND, and once using Newtonian gravity. Because 
S and the potential in the plane are the same, the Newtonian disk is supplemented with 
an inert spherical halo that gives, together with the disk, a Newtonian potential that 
equals the MOND potential of the disk alone in the plane. (The Newtonian disks have the 
same distribution functions as their respective MOND counterparts, but their Newtonian 
Q values are higher, and r-dependent, because they do not include the /i(l + L)^/^ factor 
that appears in the MOND expression for Q.) The lower the total mass of the disk is the 



- 7- 





Fig. 3. — r/l^l as a function of the disk's mass for MOND models with constant Q = 
1.5, 2.0, 2.5. The horizontal line at T/IT^I = 0.31, marks the value used in our models 
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stronger is the departure from Newtonian gravity, resulting in an increase in the relative 
contribution of the halo. In Figure ^ are shown the relative contributions of the disk to the 
total radial force as a function of radius for the Newtonian-plus-dark-matter cases. 

In order to make a quantitative comparison between the growth rates of the unstable 
modes of the different mass models we scale the time step in the simulation in proportion 
to a natural dynamical time of the model. In Figure |^ are plotted, for each model, the ratio 
of the angular frequency of the M = 0.005 mass model to to that of the specific model. As 
can be seen from the graph this ratio depends somewhat on R. We have chosen to scale the 
time step in proportion to the orbital time at i? = 0.5 (scaling the time step in proportion 
to orbital time at i? = 0.25 does not change the results qualitatively). 

The development of the instability is traced in the time dependence of the fraction of 
the disk's mass in the m=2 Fourier component of the surface density. This turns out to 
have a period of exponential growth. We take the exponential growth rate as a measure of 
the instability's strength. In Figure ^ we plot the growth rates as functions of mass for both 
the MOND and the Newtonian-plus-DM models. These are also given in Table |^ together 
with the Q value and the halo mass (for the Newtonian counterparts) of the different 
models. The growth rates are calculated using the scaled time units i.e. the real growth 
rate equals the growth rate that appears in the graph times fim(-R = 0.5)/fio.oo5(-R = 0.5) 
as was described previously. 

4. Conclusions 

From the results presented in Figure ^ we see that exponential disks having a given 
fraction of their kinetic energy in random motion, and a constant Q profile are locally more 
stable in MOND than in Newtonian dynamics, as reflected in the fact that for the same 
value for T/|W^| the MOND disks have a higher Q value. This is in agreement with the 
general result ( [Milgrom 1989| ) regarding the local stability of disks in MOND. One can see 
that the change in the dynamics occurs when one crosses over from the Newtonian regime 
to the MOND regime i.e. when the acceleration in the disk become of the order of oq. The 
added degree of stability is limited (the change in Q saturates deep in the MOND regime). 
This stems from the fact that at both the Newtonian limit and the deep MOND limit the 
equations governing the evolution of the system obey simple (but different) scaling laws. 
The basic physical mechanism behind the added stability is the relatively weaker response 
in the potential to a given perturbation in the surface density when one is in the MOND 
regime. Roughly speaking in MOND oc p and therefore 6a/ a = 6p/2p, while in the 
Newtonian case Sa/a = Sp/p, and we see that approximately a factor of two is gained in 
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Fig. 5. — The ratio of the angular frequency of the m = 0.005 mass model to the angular frequencies 
of all the mass models as a function of radius 



- 10 - 



stability. 

From the results presented in Figure ^ and Table |I] we see that the global stability of 
the disk behaves in a way similar to the local stability. As one moves from the Newtonian 
regime to the MOND regime the growth rate of the m = 2 mode (in dynamical-frequency 
units) decreases. At first (down to M ~ 0.2) the effect of MOND is similar to that of the 
added inertial halo. Below that the degree of stability continues to increase, but not as 
fast as that of the Newtonian disk-plus-halo; and it saturates in the deep MOND limit. 
In contrast, the Newtonian disk-in-halo becomes increasingly stable in the limit. The 
saturated global stability given the disk by MOND, is similar to that given a Newtonian 
disk by an inert halo with a mass that is 2-3 times the mass of the disk up to i? = 5 scale 
lengths. These results support the idea that pure MOND disks with high surface densities 
are less stable than those with a lower surface density both globally, and locally. This 
provides a possible explanation of the Freeman law as discussed in the introduction. It 
must, however, be appreciated that we cannot be sure that actual LSB galaxies are more 
stable than HSB ones because we do not know that they all have similar T/W values, as 
used in our comparison. 

Our aim in the paper has been to compare the stability properties in MOND of disks 
with different accelerations. In this, the Newtonian models have served as references so that 
the added MOND stability could be described in terms of an added inert halo. But, what 
is the significance of the comparison between the MOND and Newtonian disks as regards 
true galaxies? The disk-plus-inert-halo models we use are not what MOND predicts for 
a galaxy. If MOND is correct than in the low acceleration limit there should also appear 
to be much disk dark matter. In the present paper we have ignored the structure and 
motion in the z direction, perpendicular to the disk. A Newtonian model that will give 
the same 3-dimensional disk distribution function as a MOND pure disk, would have much 
disk dark matter that is not inert, but responds to disk perturbations. Put differently, 
MOND predicts that the dynamically determined surface density of LSB galaxies will 
be much higher than the observed surface density. When this surface density is used in 
calculating the Newtonian Q value, as it should, a much lower Q value will result in general 
( PVIilgrom 1989|) . Inasmuch as we neglect the z-structure and take the disk as infinitely 



thin the exact value of the Newtonian surface density that gives the same potential field 
as MOND is Sat = S//!"*", where /i"*" is the value of fi for the local acceleration just above 
the disk (because at every point on the disk, and at all times during its evolution, we have 
fi~^dn(f) = 27cGT,, while 2ttGTii\i = (9„0). The net result is that even in the deep-MOND limit 
MOND disks are expected to be somewhat more stable than the Newtonian disks that have 
the same distribution function (and which thus have the same r- and z-structure) because 
of the (1 + L)^/^ factor. We do not include here simulations for such Newtonian models. 
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A. Model construction 

The problem of finding a distribution function (DF) f{E, L^) can be made well posed 
for numerical solution by formulating it as a constrained optimization problem (see [Binney| 



and Tremaine 198^ ). One wants a DF that satisfies the following physical requirements as 



constraints 

/>0, 

//rfv = S(r), (Al) 
/Krfv///rfv = a2,(r), 

and, as an additional auxiliary constraint that will assure uniqueness, maximizes a certain 
functional of / such as the Boltzmann entropy or some measure of smoothness. A very 
similar approach is to minimize a single functional of / that is the sum of the errors in the 
surface density, the radial velocity dispersion, and the entropy, subject to the constrain that 
/ > 0. We have used the latter approach. We take the disk to be a finite disk whose surface 
density vanishes for r > 1. The three input functions: the potential, the surface density, 
and the desired radial velocity dispersion, are represented by one dimensional arrays (pi, Sj, 
and ayn, respectively, at the equidistant grid point rj = i/N {0 < i < N). 

Using the variables 

X = ^-E + <f>{l), Y=^, (A2) 

and given a distribution function /, the surface density and radial velocity dispersion runs 
take the forms: 

S'(r) = 1 1 f{X,Y)W{X,Y,r)dXdY, (A3) 

a'^lir) = ^'-\r)J J f{X,Y)W{X,Y,r)Z{X,Y,r)dXdY, (A4) 

W{X,Y,r) = Y-^/^r^[(l){l) - (l){r) - X] - {1 - r^)Y}-'/^ 
Z{X,Y,r) = v^,{r) = 2[Y{l-r-^)-X + (f){l)-(f){r)] 
0< X < 0(1) - 4>{r) 

2 

0< Y <-i_^[0(l)-0(r)-X]. 



Note that we work here with the choice R^ut = 1; more generally we would have to replace 
Lz in the above expressions by Lz/Rcut- 



Before discretizing the equations we make a change of variables from X, Y to r^m, ^max, 
which denote, respectively, the pericenter, and apocenter of an orbit with given energy 
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and angular momentum. The transformation from r mini max coordinates to the X,Y 
coordinates can be obtained by solving the two equations: 



^2 



Y ^ (A5) 



mm 

2 



Y = ^ma. [0(1)_0(^^^^)_X]. (A6) 

1 — rf: 



max 



From the expressions X [r mini ^ max) and Y(rminii^max) "we calculate numerically the Jacobian 
dix,Y) — ^ g^j^^ rewrite the integrals (|A3|, |A3) using the new coordinates where now the 

d{rmin,rm.ax) ^ 1' <- 3/ to 

limits of integration are: 

< Tmin < r, (A7) 
r < rmax<l. (A8) 



We also replace S'(r) in equation ( [A^ ) by the desired surface density, since the two 

become identical when a solution is found. We discretize the DF on a Cartesian grid where 

fjk f max j 1 miuf;') (-'^9) 
Tmax, = {]-l)/N 
TmiUk = {k- l)/N 

and the value of f in between grid points is defined using bilinear interpolation. Since 
interpolation is linear in fjk and the integrals in eqs. ( |A3| , |A4| ) , after the replacement of S', 



are also linear in fj^i we can obtain expressions of the form: 



a: 



n = J2^^JkfJk1 (AlO) 

■ ^Bijkfjk/^ii (All) 



/2 



where Aijk and Bijk are calculated by numerically integrating 

d(x,Y) — -W(X,Y,r) and 9(x,y) — \y{x,Y,r)Z(X,Y,r) respectively, over the relevant 
volume for each grid point. We are left with the discrete problem of minimizing the 
expression 

N ^2 

ai E (S. - + «2 E « - O + «3^[/] (A12) 

i=l i=l 

with respect to the variables fjk given Sj and a^^,. The functionals of / that we have 
used are the Boltzmann entropy which is defined as S* = — / /log/rfxdv, or a measure of 
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smoothness which we took as the norm of the gradient of f in the r^m, r^ax coordinates. 
After discretization we are left with an expression of the form 

S = T.Cjkfjklog{f,k) (A13) 

j,k 

for the Boltzmann entropy, and an expression of the form 

S = T. Djklifjk - kf + {fjk - /, k+if] (A14) 

for the smoothness functionaL We minimize ( |A12|) using an iterative scheme where at each 
step we make a Gauss-Seidel relaxation sweep, we sweep over the grid, and at each grid 
point we set a new value for fjk in such a way that it will minimize eq. (|A12|) using a 



quadratic approximation obtained from the first and second derivatives of eq. ([A12| ) with 
respect to fjk- After each relaxation sweep we decrease the weight 03 by multiplying it by 
a number smaller than one. In this way we let 03 tend towards zero as the calculation 
progresses. In Figure |^ we plot as solid curves the input constraints of the surface density 
and the square of the radial velocity dispersion and as points the values calculated from 
a numerical solution found for the distribution function. The galaxy model is a smoothly 
truncated exponential disk having a total mass of 0.005, a constant Q = 2.55, and obeys 
MOND. These models are described in section |^ As can be clearly seen from the graph the 
relaxation converges to an accurate solution. 

In Figure |] we plot the distribution in phase space using r^m — ^max coordinates of 
500,000 particles according to a distribution function found for the model. In Figure ^ 
we plot the distribution in phase space using rmin — fmax coordinates of 500,000 particles 
according to a distribution function found for the model. 



B. An N-body simulation and initial conditions 

The nonlinearity of gravity in MOND prevents one from using most standard potential 
solvers, at least in a straightforward manner. Since we have written a multigrid potential 
solver it is a natural choice to use the particle-mesh algorithm, as described, e.g. in [Hockneyl 



and Eastwood (1988)| , for N-body simulations. At each time step the density is interpolated 
from the particles to the grid; then we solve for the potential on the grid and interpolate the 
forces computed on the grid to the particle's location in order to integrate its equations of 
motion. We use the cloud in a cell (CIC) charge assignment and multi-linear interpolation 
for the force calculation at the particle's location; this algorithm is relatively fast. The same 
program can perform a simulation using MOND or Newtonian dynamics. 
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-Newionian + Halo 
-MOND 



Fig. 6. — The growth rate, in units of the dynamical frequency, for the m=2 mode as a function 
of the total mass of the disk. □ MOND, A Newtonian + Halo. 




Fig. 7. — The surface density and the square of the radial velocity dispersion for the m = 0.005, 
Q = 2.55, MOND, truncated exponential disk. The solid line represents the model and the dots 
are the numerical values calculated from the distribution function. 
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Fig. 8. — The distribution of particles in phase space for a Q=2.55, MOND, truncated exponential 
disk model having a total of 500,000 particles. 



The potential solvers and the N-body code were extensively tested using Newtonian 
dynamics and MOND. Once the potential solver has been tested as described at the end 
of section 2 and found accurate there is no difference between Newtonian dynamics and 
MOND in the rest of the N-body code. The N-body code was tested by running stable 
King models, both Newtonian and MOND models, and observing the stationarity of 
the different quantities such as the size of the system, average velocities, total angular 



momentum, linear momentum, energy, etc. [Kalnajs (1978)| reported the eigenfrequencies 



of the dominant bisymmetric eigenmodes of the isochrone/mfc models. [Earn and Sellwod| 



1995)1 used Kalnajs's distribution functions for the isochrone/12,9 models to compare the 



results they got from their expansion code to the analytic results of Kalnajs. We have run 
a simulation using our code, Newtonian dynamics, and their initial conditions. We then 
performed the same fit as they did for the pattern speeds and growth rates of the unstable 
modes. The fit between the numerical results and the analytical results were good (about 
10-15% accuracy). 

The importance of a careful initial setup for an equilibrium model is well documented 
( ^ellwood 1983| , [Sellwood 1987|) . There are two separate aspects to this: suppression of 



particle noise and choosing coordinates from the desired distribution function. Initial 
positions picked randomly produce shot-noise density fluctuations on all scales. For initial, 
near-equilibrium models the initial behavior is dominated by the collective response to the 
artiflcial noise, and this can mask the dominant modes of the continuous system ( |Sellwood 
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1983), in which we are interested. Such initial noise can be suppressed by arranging 
the particles regularly. This results in a discrete noise spectrum with large amplitudes 
at the wavelength of the particle spacing, which must be suppressed during the force 
determination. This would give a particle distribution that behaves as a smooth fluid. To 
this end, for our polar grid we place particles on rings spaced in radius according to the 
required surface density law. The number of particles on each ring must be related to the 
number of azimuthal Fourier harmonics rrimax that enters to the force. To prevent coupling 
of modes through aliases, 2{mmax + 1) particles are needed on each ring. Quiet starts are 
also possible for warm stellar disks, but it is not practicable to suppress both radial and 
azimuthal density variations at the same time ( |Sellwood 1983| , [Sellwood and Athanassoula 



1986| ). Noise in the azimuthal forces must be suppressed since we are interested in 
non-axisymmetric instabilities. One then gives the particles on the ring identical velocity 
components so that the initial orbits remain congruent. 

Choosing integrals of motion for each particle at random from the distribution function 
will result in statistical fluctuation about the intended function. These can be eliminated 
by choosing integrals for each particle in a deterministic manner such that their distribution 
is as close as possible to the required form. For example, one can use energy and angular 
momentum as the independent variables ( [Sellwood and Athanassoula 19861) ; however, any 



other set of isolating integrals in which the distribution function can be expressed would 
work equally well. 

Since we are interested in making a quantitative and systematic comparison between 
the stability of bare disks obeying MOND, and Newtonian disks with dark halos we want 
to minimize the statistical noise and employ a quiet start technique. As discussed above, 
one needs to use only a selected number of azimuthal Fourier components of S in the force 
determination. In Newtonian dynamics this is justified for linear stability analysis since 
the Poisson equation and the linearized collisionless Boltzmann equation do not couple 
modes with different azimuthal frequencies. The MOND field equation can be linearized 



around the solution of the unperturbed axisymmetric disk (as discussed in Milgrom 1989 ) 
and together with the linearized collisionless Boltzmann equation have the property that 
unstable modes with different azimuthal Fourier components are uncoupled. Instead of 
using the linearized MOND equation we use the full MOND equation, but leave only the 
desired Fourier components in the surface density that is assigned to the grid. In setting 
up the initial conditions we use the following procedure: We take the numerical solution 
obtained for the distribution function and interpolate it to a finer grid. We then calculate 
the number of particles that should reside in each cell given the total number of particles. 
This number is usually not an integer, we take the integer part and distribute the particles 
uniformly in the cell. The remaining fraction is interpreted as the probability for an 
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additional particle to reside in this cell. We then draw cells at random according to their 
relative probabilities and place at most one additional particle in a cell. At the end of this 
stage we have a list of the {rmin, ^max) coordinates of the particles. We now need to assign 
the phase-space coordinates for each particle i.e. r, 9, Vr, vg. (Here we use a polar grid in the 
mid-plane as an auxiliary for computing the discretized density distribution that serves as 
input for the MOND field equation.) We draw randomly the radial position of the particle 
taking the probability density of finding the particle at radius r as being proportional to 
v^^. If we decide to use the Fourier filtering we draw at random the angle 9 and place 
2{mmax + 1) particles at angular spacings of T^immax + adding a small random angular 
shift to each particle to seed the unstable modes. If we do not use the Fourier filtering we 
place the particle at a random angle 9 requiring that the surface density produced on the 
grid that is used by the potential solver will be as smooth as possible. 
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m 


Q 


time step 


Growth rate 


halo mass 






scaling 


MOND 


Newt+DM 


at R=l 


0.005 


2.55 


1 








0.01 


2.5 


0.84 


0.4 






0.02 


2.4 


0.7 


0.43 






0.04 


2.25 


0.58 


0.46 


0.09 


0.18 


0.08 


2.0 


0.48 


0.51 


0.36 


0.23 


0.16 


1.79 


0.39 


0.62 


0.53 


0.28 


0.32 


1.62 


0.3 


0.8 


0.8 


0.31 


0.64 


1.53 


0.22 


0.94 


0.94 


0.31 


1.28 


1.5 


0.16 


1.0 


0.97 


0.27 



Table 1: The growth rate, in units of dynamical frequency, for the m — 2 mode, and model 
parameters for the different mass models. 



